Stochastic Model for the Interaction of Buckling and Fracture 
in Thin Tension-Loaded Sheets 
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We introduce a model of fracture which includes the out-of-plane degrees of freedom necessary to 
describe buckling in a thin-sheet material. The model is a regular square lattice of elastic beams, 
rigidly connected at the nodes so as to preserve rotational invariance. Fracture is initiated by dis- 
placement control, applying a uniaxial force couple at the top and bottom rows of the lattice in 
mode-I type loading. The approach lends itself naturally to the inclusion of disorder and enables a 
wide variety of fracture behaviours to be studied, ranging from systems with a simple geometrical 
discontinuity to more complex crack geometries and random cracking. Breakdown can be initiated 
from a pre-cracked sheet or from an intact sheet where the first damage appears at random, and 
buckling sets in when a displacement vector containing out-of-place components becomes energeti- 
cally favourable over one which does not. In this paper we only consider center-cracked sheets with 
no disorder and include some results relevant to the force- and displacement-fields, and the buck- 
ling response ratio. Rather than carry out a comprehensive study of such systems, the emphasis 
presently is on the development of the model itself. 

PACS numbers: 81.40.Jj, 62.20.-x, 05.40.-a 



I. INTRODUCTION 

Understanding how, and when, materials break are im- 
portant in many engineering applications. This is so for 
a number of reasons - the motivation to study fracture 
may, for instance, be related to safety issues, such as 
determining when cracks form in concrete structures, or 
it may be one of economical gain, as in the case when 
the runnability of a printing press in the paper industry 
is considered. Few materials, be they natural or man- 
ufactured, are perfect, however. Hence, the disorder in 
the micro-structure needs to be accounted for in order to 
obtain a realistic description. 

Over the past fifteen years methods have emerged 
within the statistical physics community to successfully 
tackle just such problems pj. These methods are in- 
termediate between the microscopic, or first principles, 
approach and the mean-field type of approach. In the 
former case fracture properties are derived from inter- 
molecular or inter-atomic forces, representing a problem 
which is both theoretically demanding and heavy on nu- 
merical resources. In the latter case disorder cannot be 
included in a satisfactory way. This is a big drawback 
since the presence of disorder in a material is crucial to 
the way it fractures. Disorder affects the stress field in 
such a way as to enhance the already existing hetero- 
geneities. This interplay, between a constantly evolving 
non-uniform stress field and local variations in material 
properties, can nevertheless be handled in a numerically 
tractable way using lattice models. 

The most common lattice models used in engineer- 
ing applications are finite element methods (FEM), the 
implementation of which is usually based on commer- 



cially available computer codes. The lattice models cur- 
rently used in statistical physics differ somewhat from 
the FEM-approach in that the grid used is regular, i.e., 
the same everywhere, rather than one which adjusts the 
mesh size according to where the stress field is most in- 
tense. Although FEM modeling is certainly more suit- 
able in describing homogeneous materials, the require- 
ment that the stress field should vary slowly over each 
element makes the approach cumbersome in the pres- 
ence of heterogeneities. In the stochastic lattice model, 
however, the nodes are thought of as being connected 
by objects such as elastic beams or current carrying el- 
ements. While in some respects being less sophisticated 
than FEM methods, the interpretation of the algorithm 
is much more transparent and the approach also has the 
advantage of allowing disorder to be included quite gen- 
erally. 

In the stochastic models, the local equilibrium in force 
and moment is considered on a mesoscopic scale, i.e., on 
a scale much smaller than the external dimension of the 
lattice but still sufficiently large for the forces to be gov- 
erned by well known physical laws. In this sense it is 
also a very good alternative to the far more complicated 
approach of including disorder on the microscopic level. 
Since only the nearest neighbours on the lattice are in- 
cluded the calculation of the displacement field reduces to 
the inversion a sparse matrix, enabling reasonably large 
systems to be handled computationally. 

One feature which is of a phenomenological nature, 
however, is the breaking rule - the choice here is guided 
by intuition rather than by the inner workings of the 
model itself. In other words, breaking does not arise as a 
natural consequence of the calculations. This is actually 
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an advantage in that the mechanism by which the sys- 
tem ruptures can be tailored to suit different engineering 
requirements. If we regard thin planar materials, for in- 
stance, the energy required to propagate a crack across a 
given area is usually much lower in mode-Ill fracture, i.e., 
tearing, than in the pure tensile loading of mode-I frac- 
ture. Familiar examples of disordered materials which 
behave this way are textiles and paper. 

Most of the research done so far in stochastic lattice 
modeling aims to identify the underlying general princi- 
ples of the fracture process rather than to address tra- 
ditional problems in fracture mechanics. In this paper 
the plane beam model Q, which has been used previ- 
ously to study scaling laws in fracture, is extended to 
include a specific, and practical, aspect of fracture which 
is very important for thin sheet materials, i.e., buckling. 
As is well known, buckling can profoundly influence the 
residual strength of such materials 0, Q El ■ But before 
devoting our attention to this problem in full, we briefly 
mention part of the background which has inspired the 
use of lattice modeling as a tool in statistical physics. 

In modeling experiments of random media, the fea- 
ture which by far has received the most attention is 
the morphology of crack surfaces. Many surfaces in na- 
ture are found to be self-affine, i.e., statistically invari- 
ant with respect to anisotropic scale transformations. 
The morphology of such surfaces can be described by 
simple scaling laws, behaving very much like fractal ob- 
jects [6|. These scaling laws provide a theoretical frame- 
work whereby much information can be summed up in 
a few parameters. Certain features have been found to 
share a common basis with other, seemingly unrelated, 
problems such as deposition and growth processes, or 
transport properties in random media 7] . In the case of 
fracture it has been established that crack surfaces scale 
as W ~ , where L is the system size, W is the rough- 
ness and C is the roughness exponent. Other scaling laws 
have been studied, e.g., in connection with the distribu- 
tion of stresses, or for the total amount of damage found 
at various stages in the breakdown process. 

By far the most popular tool in such studies has been 
the random fuse model |8j . In the fuse model, the nodes 
on the lattice are connected by current-carrying elements, 
i.e., fuses. The threshold for the amount of current which 
may flow through each fuse is chosen from a random dis- 
tribution. Hence, in the breakdown process a fuse is ir- 
reversibly removed from the lattice once its threshold is 
exceeded. A new distribution of currents is then calcu- 
lated before the next fuse is removed, and so on, until an 
uninterrupted path can be traced across the system. Al- 
though it really describes electrical breakdown, the fuse 
model is often referred to as a scalar model of fracture due 
to the similarity in form between Ohm's law and Hooke's 
law of linear elasticity. Results obtained for £ with the 
fuse model are found to be different in two and three di- 
mensions, however. The results are £ = 0.74(2) in two di- 
mensions 9] and £ = 0.62(5) in three dimensions Al- 
though the former seems to agree with experimental find- 



ings, the latter does not. Furthermore, the type of forces 
involved on the meso-scale also seem to make a differ- 
ence, i.e., the results obtained with a scalar model differ 
from those obtained with a vectorial model. Specifically, 
in calculations with the elastic beam model £ = 0.86(3) is 
obtained in two dimensions The difference between 
the results of the two and three dimensional fuse model 
indicates that the additional degrees of freedom afforded 
by the (three-dimensional) buckling beam model should 
provide a lower estimate for the roughness exponent in 
the vectorial problem as well. Since the observed value 
in real materials, i.e., C = 0-8 m f ac t does lie below 
the two dimensional beam lattice result, it would be in- 
teresting to see if the buckling beam lattice reproduces 
the universal value observed in nature. 

However, although such fundamental aspects of the 
fracture process are certainly interesting, the subject of 
how buckling affects the scaling laws are left for future 
study. The focus in this paper is instead on the devel- 
opment of a lattice model that realistically includes the 
buckling behaviour observed in thin sheet materials. The 
characteristic out-of-plane deflection known as buckling 
is perhaps most frequently associated with thin plates 
or beams under compressive loading. Presently, how- 
ever, we concern ourselves with the special case of a thin 
planar structure under tensile loading. The interaction 
of buckling with fracture in such circumstances is a well 
known phenomenon, although it has often been neglected 
in fracture mechanics analyses due to the extra compli- 
cations involved. One of the characteristic features of 
buckling in a thin tension-loaded sheet is that a stable 
out-of-plane configuration is obtained after buckling has 
set in. This is in stark contrast with the case of com- 
pressive loading, where loss of stability usually signals 
complete breakdown. 

Practically all previous work considers the effect buck- 
ling has on the strength properties of an already cracked 
plate or a plate with a geometrical discontinuity such as 
a circular hole or a rectangular cut-out. If the physi- 
cal parameters of the plate are such that buckling can 
be expected before the crack begins to grow, the residual 
strength of the plate will be significantly lower than what 
would otherwise be expected, based on an analysis which 
does not take account of buckling. The present study of 
fracture and buckling will be more general in scope. In 
other words, we also regard sheets which, in their initial 
state, have no cracks or other discontinuities. Instead, 
cracks form by a complex process which depends on the 
evolving distribution of stresses and its interaction with a 
disordered meso-structure. The onset of buckling in this 
scenario, and the effect buckling has on the fracture prop- 
erties, will vary according to the type of disorder used, 
i.e., weak or strong. Nonetheless, even for weak disorders 
the final crack which breaks the system will only rarely 
appear at the exact center of the sheet, and even then 
the situation will usually be complicated by additional 
cracks in the vicinity - cracks which interact with the 
main crack so as to alter the distribution of stresses and 
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FIG. 1: On the left-hand side is shown a lattice of size L = 5 
where a force couple has been applied uniformly on opposite 
edges. The strain imposed is consistent with mode-I type 
fracture and corresponds to a displacement 5L — 1 in the 
F-direction. The enumeration scheme of the neighbouring 
beams is shown on the right-hand side, where a rotation at 
node i (center dot) induces shearing forces and bending mo- 
ments in the neighbouring beams. 



hence also the exact shape or mode of buckling. 

The emphasis here, however, is on the development of 
the model itself. For illustration purposes, a few results 
are included on uniform systems with a center-crack. In 
section ITT1 the plane beam model is briefly reviewed, be- 
fore the equations describing the out-of-plane behaviour 
are derived in section IIIII Typical stress and displace- 
ment fields are shown in sections ll V I and IV! respectively, 
before the initialization of buckling is discussed in sec- 
tion IVII and a fracture criterion defined in section IVIII 
where results for the buckling response ratio of a center- 
cracked sheet are included. 



II. PLANE BEAM LATTICE 

The beam model may be defined as a regular square 
lattice of size L x L, where the spacing is one unit length, 
and each node in the horizontal and vertical in-plane di- 
rections is connected to its nearest neighbours by elastic 
beams. A beam is then fastened to other beams in such 
a way that, upon subsequent displacement of neighbour- 
ing nodes, the angle between beams remains the same 
as in the original underlying square lattice, see Fig. ^ 
Furthermore, all beams are imagined as having a certain 
thickness, providing finite shear elasticity. 

Beginning with the simple two dimensional beam 
model, there are three possible degrees of freedom, i.e., 
translations in the horizontal (a;) and vertical (y) direc- 
tions, and rotations about the axis perpendicular to the 
plane (w). As shown in Fig.^ this allows for both bend- 
ing moments and transverse shearing forces, in addition 
to the axially tensile, or compressive, forces. 

For any node i, the nearest neighbours j are numbered 
in an anti-clockwise manner, beginning with j = 1 to the 
right of i. Defining Sr — Vj — Ti, where r 6 {x, y, u>}, the 



forces on i due to j = 1 are 
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for the moment due to angular displacements w, shear 
and transverse force due to displacements y, and axial 
strain due to displacements x, respectively. Expressions 
for j > 1 are analogous. 

Prefactors characteristic of the material and its dimen- 



sions are 
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where E is Young's modulus, p and I the area of the beam 
section and its moment of inertia about the centroidal 
axis, respectively, and G the shear modulus [l3T |. 



The conjugate gradient method 
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are the components of force and moment. 

The material is assumed to be brittle, i.e., each beam 
is linearly elastic up to the breaking threshold. Using 
tA and tM for the thresholds in axial strain and bending 
moment, respectively, a good breaking criterion, inspired 
from Tresca's formula, is 



-V 



\M\ 

tM 



> 1, 



(9) 



where \M\ — max(|Mi|, \Mj\) is the largest of the mo- 
ments at the two beam ends i and j. 

The fracture process is initiated by imposing an ex- 
ternal vertical displacement which at the top row corre- 
sponds to one beam-length, i.e., 6L = 1, see Fig.^ The 
lattice now consists of horizontally undeformed beams 
and beams which in the vertical direction are stretched 
lengthwise. The first beam to break is that for which the 
ratio A/t a is largest, this being the vertically oriented 
beam which has the lowest value of tA- If ah thresh- 
old values are the same, the next beam to break will 
be one of the nearest lateral neighbours since these now 
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FIG. 2: A disordered beam lattice of size L = 19, which is 
strained to failure in mode-I type fracture, i.e., by applying 
a force couple at the top and bottom edges. The presence of 
a central crack leads to the build-up of compressive stresses 
around the crack edges, causing the structure to deflect out 
of the initial rest plane. 



carry a larger load than other beams on the lattice. The 
case of no disorder is thus one in which the crack prop- 
agates horizontally from the initial damage, taking the 
shortest possible path to break the lattice apart. In- 
troducing disorder in the breaking thresholds, material 
strength varies across the lattice and consequently the 
crack will not necessarily develop from the initial dam- 
age point. Instead microcracks and voids form wherever 
the stress concentration most exceeds the local strength, 
i.e., wherever Eq. (JjJ dictates that the next beam should 
be broken. Towards the end of the breakdown process 
smaller cracks merge into a macroscopic crack, forming 
a sinuous path which ultimately traverses the width of 
the lattice and thus breaks it apart, see, for instance, 
Fig. El In this scenario the quenched disorder on the 
thresholds and the non-uniform stress distribution com- 
bine to determine where the next break will occur. The 
stress distribution itself also continually changes as the 
damage spreads. 

Throughout the process, the equilibrium stress field 
is re-calculated by use of Eq. J^J each time a beam is 
removed. The stress field therefore relaxes at a rate much 
faster than the process by which the crack grows. Hence, 
the model describes quasi-static fracture. 



III. BUCKLING BEAM LATTICE 

The displacements of a real material, even if its ge- 
ometry is essentially confined to a plane, will generally 
occupy three dimensions. For instance, when opposite 
forces are applied uniformly along the top and bottom 
edges of a sheet of paper, with the object of straining it to 
failure, significant displacements will be observed in the 
direction perpendicular to the sheet, see Fig. [5] This be- 
comes especially evident wherever sizable cracks appear. 
Reasons for this behaviour are deviations in the sym- 
metry of the material itself, or its properties, about the 
plane through which the externally applied forces act. In 
some cases such deviations may simply correspond to an 
uneven thickness, or they may be caused by local varia- 
tions in density, a gradient in the orientation of the micro 



structure, and so forth. 

To include this behaviour, the plane beam model must 
incorporate at least two additional features. One is the 
random variation of the material in the out-of-plane di- 
rection. Since lattice modeling reduces the material to a 
set of points corresponding to the nodes on a mathemati- 
cally precise two-dimensional lattice, the most convenient 
approach is to impose a very small randomly chosen verti- 
cal displacement on each node. This is discussed in more 
detail in section IVTl The other feature to be included, 
and the topic of the present section, is the physics of 
the forces which create, and maintain, the out-of-plane 
displacement field. 

In the buckling beam model we have one translational 
and one rotational displacement relevant to each of the 
principal axes, i.e., six degrees of freedom, with the ma- 
trix system 



(10) 



replacing Eq. (J5J. Presently the forces are projected onto 
the XY-, XZ and FZ-planes, and hence Xi of Eq. i(TU)l . 
that is, 
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where the term ^XZ^ 1 ', for instance, is the x-component 
of the buckling (B) force due to j = 1, as projected onto 
the XZ -plane. Axial and transverse contributions are 
denoted (A) and (T), respectively. 

The rotational displacements about the Y- and X-axes 
are denoted u and v, respectively, and z is used for verti- 
cal displacements along the Z-axis. A coordinate system 
is placed on each node, whereupon forces and moments 
are expressed as functions of the displacements. To this 
end, an elastic beam with no end restraints 0] is consid- 
ered, as in the case of the plane model. In the buckling 
model, however, the coordinate system is additionally ro- 
tated about the relevant angle within the XZ-, Y Z- or 
XY-plane, i.e., u, v or w. 
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FIG. 3: The buckling term, —^P8, at node i due to j = 1 in 
the case of an axially compressive load, showing the angular 
displacements, Ui and Uj, the bending angle, 9, the axial force, 
F, and the component P, of F, which is parallel to the beam 
axis at node i. Also shown is the original XZ-system and the 
X'Z'-system. 



With the exception of the signs on the terms of 



Eq. i|12|) , contributions from neighbours j = 1 and j = 3 
are similar, as are those from j = 2 and j ' = 4. 
Consequently, if we define 

Pi = 5 [!-(-!)'], (13) 
with <7j = 1 — j»j , and 

j-i 

ri =n(-ir (14) 

n=0 

with Sj = (— l) J Vj, for notational convenience, then the 
total force on z along the X-axis, with the contributions 
from all four of the neighbouring beams having been in- 
cluded, reads 
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In Eq. (|15fl . moreover, 
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is an angular correction to 
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the latter being the projection onto the XZ-plane of the force along the axis of the beam. If we consider the j = 1 
component in Eq. (|12|) . 
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is the contribution due to elongation or compression along the axis of the beam, 
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is due to forces which are transverse to the axis of the beam, and 
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is the contribution due to buckling. The latter arises when a beam in a bent configuration is under compressive 
or tensile axial loading, see, e.g, Fig. where the term ^XZ^ has been shown. Out-of-plane contributions with 
components along the X-axis are 



\XZ, 
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12 
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sin m 
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from transverse forces, and 



r> — cos ou sm Ui 
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from buckling. 

The buckling term, as obtained in the lowest order ap- 
proximation, is essentially the product of a bending an- 
gle, 9, and an axial force component, P, the latter being 
parallel to the axis at the opposite end of the beam. The 
component of the buckling reaction in the X'Z'-system 
which lies along the X-axis in the XZ-system is then 
^XZf\ i.e., Eq. C3. 

Eq. (|2ip. moreover, corresponds to that part of the 
transverse force (including shear) which does not include 
buckling and is similarly obtained, i.e., by rotating the 
axes in Eq. J2J. 

Finally, the axial term becomes 
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when the out-of-plane displacements are set to zero. 
Hence, in this case, only when the rotation of the XY- 
system onto the X'Y'-system is neglected does Eq. lfT^|) 
reduce to Eq. ©, of the plane beam model. 

Although forces and displacements on a beam under 
simultaneous axial and transverse loading cannot, in gen- 
eral, be obtained by superposition, combinations such as 



1XZ\ '+^XZr } in Eq. (JT2J result when only the leading 
terms, in P, are retained after inverting the expressions 
of Ref. . This also causes the buckling term in tensile 
loading to be the same as that in compressive loading, a 



change of sign being the only difference. 

The expression for the Y"-component is similar to that 
of the X-component, and is obtained by changing around 
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The Z-componcnt, furthermore, is 
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i.e., also similar in form to Eq. I|12|) but with lines number 
two and five omitted. 

The full expressions are then 
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for the force in the horizontal direction, and 
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for the force in the direction perpendicular to the rest plane. 

Considering next the rotational contributions, a beam under axial loading, which is simultaneously bent, is shown 
in Fig. 0] In this case a buckling term arises which is again the product of a bending angle, 9, and an axial force 
component. For rotations about the Z-axis, this gives 



Wi = r^ (1) + lXY^ + ri^ (2) + lXY^ + «ir< 3) + lXY^ + ™XY^ + IXY^. (28) 




FIG. 4: The contribution \P9 to the in-plane moment at node 
i from j — 2, due to buckling, in the case of a tensile axial 
load. Shown are the angular displacements, Wi and Wj, the 
bending angle, 8, the axial force, F, and the component P, of 
F, which is parallel to the axis of the beam at the opposite 
end, i.e., at node j = 2. 



For each beam in Eq. H28fl there are two terms, one anal- 
ogous to Eq. and denoted (M), and one extra term, 
such as ^XYi , which is the contribution due to the 
beam being simultaneously bent while under axial load- 
ing. Similarly, we have 
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for rotations about the K-axis, and 
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for rotations about the X-axis, where (Q) denotes the 
torque. In Eq. 1|29[) . the torque is simply 

Qyy/ 2) = #u, (3i) 

where, assuming w > t, the material constant is 



(32) 



when w denotes the width of the beam cross section and 
t its thickness. The buckling term reads 
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and that part of the bending moment which does not 
involve buckling becomes 
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when the axes are rotated. Eq. I|29|l . when written out in full, is now 



1 + Su) sin Ui — Sz cos Ui + — 



(34) 



Ui 



1 



ft+^ p AiY 



Su 



l r 



(l + Tj8x) sinui — rjSz cosui + 



Su 



E 

3=1 



F^'pj— cos Su - q^Su 



(35) 
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and Eq. is analogous, i.e., 

P + 12 j =1 I ix. J J J 3 - =1 

Finally, for rotations within the XY-plane, Eq. I|28|) becomes 



Sv 

qj — cos Sv — Pj^Sv 



(36) 



4 r 



Wi = 



l—Y 

209 +5) ft 



sin to, — Sj (pjSx + qjSy^J sin lOj — rj (qjSx + PjSy) cos iOj + 



<5w 







4 1 4 

tit Z Jw - 1 Z F ^ )(5w cos 



(37) 




FIG. 5: Node i and its nearest neighbours j = 1-4, shown 
when the lattice is in an advanced state of buckling. The 
plane passing through i is uniquely defined by any j and j ± 1 
neighbours of i, as shown by the broken lines, and is no longer 
parallel to the XY-plane. The out-of-plane reaction (a) is 
normal to the X'Y'-plane and (b) is a bending moment about 
the Y'-axis. 

In the six components of Eq. (|10|) . derived above, pref- 
actors characteristic of the beam and its dimensions vary 
according to the principal axis of bending. Hence, in 
Eq. we use 

h = ^wH (38) 

for bending within the XY-plane, and 

Ix = -^wt 3 = Iy (39) 

for bending within the YZ- and XZ-planes. We have 
then assumed beams with a rectangular cross-section, as 
already noted in connection with Eq. (|32|l . This is con- 
venient in the study of how thin sheets behave during 
fracture, since one may then simply visualize beams with 
a flat profile, see Fig. El In the present calculations the 
chosen width-to-thickness ratio is 10:1, so that resistance 
towards bending within the plane is much larger than 
that which governs out of plane bending. 

In the following, results are displayed for non- 
disordered systems with a central crack. To illustrate 
the nature of the forces, moments and displacements in- 
volved, sections of the lattice parallel with the crack are 
referred to as J = 1, 2, L + 2. Hence, on the bot- 
tom part of the lattice in Fig. [3 the set of nodes i = 1 



to i = L + 1, located on the same row parallel with the 
X-axis, is referred to as J = 1. With a total of L + 2 
rows J parallel with the Y-axis, the "near" edge of the 
crack coincides with J = L/2 + 1 while the "far" edge 
coincides with J = L/2 + 1. Likewise, the set of nodes 
i = ltoi = L + 2 parallel with the Y-axis is referred to 
as, from left to right, 1 = 1 to I = L + 1. 



IV. DISPLACEMENTS 

The equations governing force and moment in a buck- 
ling beam lattice were derived in the previous section. At 
present we have not taken account of the Poisson con- 
traction which is observed in elastic systems - at least 
not at the level of the individual beam. Such an effect 
does show up, however, on length scales spanning sev- 
eral beams. Of course, in the macroscopic behaviour of 
the lattice, an example of this is precisely the buckling 
behaviour we intend to study. The bulging of the crack 
edges shown in Fig. |3 for instance, comes about as a re- 
sult of transverse compressive stresses which develop in 
the neighbourhood of the crack. 

Fig. shows the in-plane displacements Xi and yi in a 
lattice of size L = 100 at various stages of crack advance- 
ment. The out-of-plane deflection Zi(N) is also shown 
on the right-hand side. Displacements refer to the ini- 
tial coordinate system on each node. In the specific ex- 
ample shown the crack grows towards the left-hand side 
of the lattice, with the out-of-plane deflection increas- 
ing with the extent of the crack opening. As previously 
mentioned, fracture is initialized by displacing the top 
row a unit distance. In the absence of geometrical dis- 
continuities, each horizontal row J is then incrementally 
displaced by an amount (L + 1) _1 with respect to the 
previous row J — 1. In the absence of cracks, the dis- 
placement field yi(J,N) then consists of a set of equidis- 
tant lines between zero and one. With a crack present, 
this is altered into the pattern shown in Fig. HO e.g., for 
yi(0). As expected, transverse displacements Xi(N) are 
largest close to the face of the crack. If we consider 2^(0), 
and move from left to right along the edge of the crack, 
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FIG. 6: Displacement fields across the width of a lattice of size L = 100 (J = 1 to I = 101) with an initial center-crack. On 
the left-hand side Xi(N) is shown for J = 1, 3, 5, etc., up to and including the crack interface, i.e., J = 51. At center is shown 
Ui(N) for J = 1, 3, 102. On the right-hand side, Zi(N) is shown for J = 1, 3, 51, including also the far side of the crack 
interface, i.e., J = 52. The number of beams broken is N, and crack extent in the four stages shown is I — 34 — 68 (N = 0), 
I = 16 - 70 [N = 20), I = 6 - 70 (N = 30), and I = 1 - 75 (N = 40). 



beams are seen to be stretched wherever the slope is posi- 
tive and compressed wherever it is negative. As the crack 
grows the net effect, however, is to cause the lattice to 
contract in the transverse direction, e.g., with the edge 
on the left-hand side moving inward by about 1.5 beam 
lengths in the case of Xi(40). 

The rotation of axes mentioned above is necessary to 
obtain the correct feedback between the force compo- 
nents in the system, such as mutual consistency between 
XY-forces and the Z-forces. To illustrate this, regard 
the lattice before it begins to buckle, i.e., when the stress 
field is confined to the XY-plane. When a crack grows 
beyond a certain critical size, interaction between the 
stress field and random variations in the Z-direction ini- 
tiates out-of-plane displacements which ultimately result 
in a buckled lattice. The driving forces are terms normal 
to the XY-plane, i.e., terms such as f XZ^ and f YZ\ 2 \ 



which belong to and Z^ 2 \ respectively. These terms 
are not large. In the flat lattice, for instance, the last two 
terms on the right-hand side of 

X™ = $XXP +IXYP + *XY^ (40) 

+ IxzP + ^xzP 

are identically zero while the terms ^XYP^ and ®XY^ 
are very small in comparison with the leading axial term. 
The terms ^XZ^ 1 ' and ^YZ\ 2 \ however, although being 
small in comparison with either X F^ or x Tp\ are non- 
negligible. This owes to the fact that there is no phys- 
ical obstruction in the lattice to inhibit displacements 
in the Z-direction, i.e., there is no leading ^ZZi term. 
Moreover, as can be seen from Fig. [5] when the lattice 
is in an advanced state of buckling there will be regions 
where the out-of-plane buckling reaction is inclined with 
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FIG. 7: A lattice of size L — 50, with an initial center-crack, 
shown at four different stages of fracture. The number beams 
broken are, from (a) to (d), N = 0, 10, 20, and 34, respec- 
tively. 



respect to the Z-axis, resulting in contributions of the 
type *XZ[ l ] ^ 0. These are smaller than fXZ^\ but 
are also assumed to be non-negligible. In order to include 
such terms the axes are rotated to the local deflection of 
the lattice, whereupon the components along the X-, Y - 
and Z-axes are obtained. 

An example of the effect this has is when a large crack, 
perpendicular to the force couple, opens up in the center 
of the lattice, as in Fig. |2 or Fig. Although the ini- 
tial displacement of the crack edges is normal to the rest 
plane, an in-plane component appears as the crack grows, 
i.e., the near edge of the crack is pulled slightly back along 
the negative Y-axis while the far edge is pulled forward in 
the opposite direction. This can be seen clearly in Fig. El 
where in yi(20) the rows nearest to the crack edges are 
displaced below or above the fixed values of the top and 
bottom rows. In the case of the near edge this means that 
displacements are negative, i.e., they have moved slightly 
backwards with respect to their equilibrium positions in 
the unrestrained-strained lattice. These displacements 
become more pronounced as the crack grows, as can be 
seen from yi(30) and yj(40). 

It is particularly instructive to compare the displace- 
ments of a buckling lattice with those obtained for the 
same lattice when the out-of-plane degrees of freedom are 
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FIG. 8: Comparison between the crack-edge displacements, 
obtained for the buckling lattice shown in Fig. [7| (thick lines) 
and the same lattice when the out-of-plane degrees of freedom 
are suppressed (thin lines). At the top, the extent of the initial 
center-crack is J = 18 — 34. In subsequent stages, the crack 
extent is I = 12 - 38 (N = 10), I = 2 - 38 (N = 20) and 
J = 1 — 51 (N = 34). For Xi(N) the near edge of the crack is 
shown and for yt(N) both edges are shown. 



suppressed. In Fig. [7| a lattice of size L = 50 is shown 
in four stages of crack advancement. The correspond- 
ing in-plane displacements of the crack-edges are shown 
as thick lines in Fig. 03 with thin lines representing the 
same lattice in a non-buckling fracture mode. In the lat- 
ter case, the yi displacements are seen to be confined 
between the fixed values of the top and bottom row, the 
physical structure of the lattice itself effectively acting 
as an obstruction to displacements outside this range. 
A further feature that can be noted concerns the afore- 
mentioned "Poisson" contraction. This effect is seen to 
be present in a non-buckling lattice as well, although to 
a much lesser degree. As for the angular displacements 
u>i (not shown), these are seen to be somewhat larger in 
the buckling case, except for the peak values obtained at 
the crack tips, which are more or less the same. These 
peak values increase only as the crack nears the outer 
boundaries of the lattice, where Xi is large. 

In general, the expressions derived for force and mo- 
ment from Ref. |l5j| are accurate for small displacements 
only. This assumption we simply extend to all displace- 
ments. A second reason for rotating the axes, then, is 
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FIG. 9: Transverse stresses in a lattice of size L = 100 with 
a central crack which extends from I — 24 to / = 78. Shown 
are stresses in the non-buckling (A) and buckling (B) fracture 
modes, for every row J up to and including that which coin- 
cides with the near edge of the crack, i.e., J = 51. Negative 
values correspond to compressive stresses. 

to conserve the consistency of the approximations used. 
To illustrate the point one may, for instance, regard a 
straight beam which, in its rest state, lies along the X- 
axis. The term JXZ^ then expresses the transverse 
force as a function of 6z. When Ui is large, however, and 
the coordinate axes are fixed, Sz, in addition to the trans- 
verse force, also implies some measure of axial strain in 
the beam even though this effect has already been in- 
cluded via Eq. t[L7). Hence, rotating the axes precludes 
the introduction of systematic errors, e.g., due to angu- 
lar deflections of the lattice. It also improves the quality 
of the approximations used since angular displacements 
are rendered less severe in a rotated coordinate system. 
Large angular deflections usually involve a number of 
nodes on the lattice and hence by rotating the axes we 
avoid that too many errors accumulate. 

It should also be pointed out, moreover, that with three 
rotational degrees of freedom there will in principle be 
several displacement combinations which correspond to 
a given space orientation of the beam axis. Hence the 
projection of forces into the XY- XZ- and FZ-planes 
is an approximation based on the assumption that large 
deflections about more than one axis simultaneously are 
rare. Otherwise the exact orientation of the beam would 
be history-dependent, i.e., it would depend on the se- 
quence in which the (final) angular displacements Ui, Vi 
and Wi were incremented. 



V. FORCE COMPONENTS 

The out-of-planc force components are small, but their 
collective effect has a significant impact on the stress 
field. In the presence of significant cracks there is a feed- 
back from the ^-displacements which allows the XY- 



displacements of the buckled lattice to relax with respect 
to the XF-displacements of the flat lattice. One exam- 
ple of this is the transverse compressive stress stored in 
the region in front of and behind the crack in the non- 
buckling lattice. Buckling releases this stress, as can be 
seen in Fig. where the j = 1 component of Eq. (|ll(l 
is shown in the non-buckling (A) and the buckling (B) 
cases. In (A) a region of compressive stress confined be- 
tween the crack tips is seen to extend for a distance of 
about 6-8 rows away from the crack edge. In (B) only 
a vestige of this is left, and then only in the immediate 
vicinity of the crack. The tensile stress at the crack tips 
increases slightly in the buckled configuration. 

In the non-buckling beam model, the extra non-linear 
terms which arise when the beam is simultaneously bent 
while under axial compression, or tension, are of lesser 
importance. All forces now act within the structure 
which defines the plane so that, in calculating the in- 
plane displacement field, corrections such as those due to 
Wi in Eq. (|23(l may be neglected. The out-of-plane dis- 
placement field, on the other hand, is obtained from an 
equilibrium state in force and moment between a number 
of terms which are individually small. 

For instance, the axial, transverse and buckling compo- 
nents which make up the j = 1 contribution to Eq. (|25|) 
are shown in Fig.^B Here the vertical scales on the three 
subplots have been adjusted to the relative sizes of the 
components. The magnitudes, furthermore, refer to the 
scale of Fig. [5] and thus gives an idea of the "smallness" 
of the out-of-plane force components. In agreement with 
Fig. El the forces in Fig. are seen to be most signifi- 
cant within a region nearest to the crack edge, extending 
about 6-8 rows to either side of the crack. At the on- 
set of buckling the axial and transverse terms, ^XZ^ 
and ^XZ^\ are identically zero while the buckling term, 

? XZf } shown in Fig. is non-zero. In this-situation 
the sum of contributions from j = 1 — 4 is non-zero. As 
the out-of-plane deflection increases, an equilibrium is 
approached where the sum of forces is zero. At this equi- 
librium the buckling terms, e.g., terms such as ®XZ^ 
in Fig. I1UI remain non-zero. This is also the case with 
the out-of-plane moments, shown in Fig. ^] The three 
subplots are not mutually to scale in this case, but the 
magnitudes again refer to the scale of Fig. |5J Out-of- 
plane moments are seen to be somewhat larger than the 
axial and transverse buckling forces of Fig. Hence, 
at the equilibrium, the most significant of the non-linear 
terms are those relevant to the momentum, ®XZ^ being 
about five times larger than fXZ^. 

It is instructive to see what happens when buckling 
terms such as ^XZ\ 1] and ^XZ\ 1] are removed. Shown 
in Fig. El at the onset of buckling, is the movement of 
the crack-edge as a function of time, the time-steps be- 
ing defined by the iteration procedure which locates the 
equilibrium of force and moment. Just before the point 
at which equilibrium is reached, all buckling terms are 
"switched off" whereupon the remaining forces set about 
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FIG. 12: The movement of the crack-edge as a function of 
the time-steps in the numerical iteration, shown for a lattice 
of size L — 40 at the onset of buckling and for a central crack 
between 1—14 and / = 28. Just prior to the point at which 
the equilibrium, is reached, at time 7b g = 1570, all terms 
such as ?XZi are "switched off". 



FIG. 10: Out-of-plane force components. Shown are the axial 
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and buckling "XZf contribu- 
tions to Z^ of Eq. 1251 . Contributions to zf^ are similar, 

being mirror 
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but with the contours of 
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tions to and Z[^ are smaller. Lattice parameters used 
and contours shown are the same as in Fig. 



to locate a new minimum of elastic energy. This new 
minimum, of course, is none other than the flat configu- 
ration. As would be expected, not only do terms such as 
^XZ^ 1 ' and ®XZ^ cause buckling, they also sustain it 
once it has been established. 
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FIG. 11: Out-of-plane moment components about the Y-axis. 



Shown are the torsional yYY} 2) , axial fXZ\ L >, and buckling 
uXZ^ contributions to Ui in Eq. 1291 . Lattice parameters 
used and contours shown are the same as in Fig. 



Finally some remarks on the angular correction in 
Eq. (|17f) . which is included to allow for the possibility 
that axial force may increase or decrease as a consequence 
of bending. In Eq. Ijl7|l it is assumed that the additional 
elongation due to bending can be obtained from a mul- 
tiplicative factor. This factor is based on the ratio of a 
circular arc 0] to a straight line, the former being the 
semi-circle defined by the angular difference Su at the 
end-points and the latter the line which connects these. 
On the level of the individual beam, the presence of inflec- 
tion points are neglected in this approximation. In other 
words, up-down curvatures may only occur in combina- 
tions of two or more beams in an end-to-end alignment. 
Furthermore, as can be seen from Eq. (|17fl . the effect of 
in-plane bending moments, or transverse displacements 
perpendicular to the rest axis of the beam, are neglected 
as contributions which would otherwise add to the axial 
length of a beam. 



VI. INITIALIZING THE OUT-OF-PLANE 
DEFLECTION 

An important feature to be included in the model is 
the random variation of the material in the out-of-plane 
direction, as was remarked in section ITTT1 In thin mate- 
rials such as paper, cloth, membranes and so forth, the 
most important factor influencing the behaviour during 
fracture is not the three-dimensional structure of the ma- 
terial itself. Rather it is the out-of-plane deflection of 
this structure which makes a difference. Nevertheless, 
random variations in the thickness direction provide an 
important part of the mechanism which initiates buck- 
ling. This is because such variations combine with the 
externally applied force and the emerging cracks to create 
local forces and moments which are not perfectly aligned 
within the plane. Once a buckled configuration has been 
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FIG. 13: Buckling modes for ten samples of a lattice of size L = 100, with a center-crack between I = 34 and I — 68. For the 
lower half of the lattice, contours of every other row J = 1, 3, are shown, up to and including J = 51, i.e., the near edge 
of the crack. The far edge, J = 52, is also shown. The only difference between the samples is the random variation used to 
initialize the out-of-plane deflection. 



established, however, the variation in the thickness direc- 
tion is far less important. 

Since we presently regard the out-of-plane deflection of 
a structure which has no vertical extent, buckling must 
be initiated by other means. Specifically, in modeling 
the fracture process, the equilibrium stress field is re- 
calculated by use of Eq. 11U|) after a beam has been re- 
moved. At each step of this process, i.e., for each beam 
removed, a sample-specific random noise in the form of a 
small vertical displacement is imposed on all nodes of 
the lattice. Presently, we use a random number uni- 
formly distributed on the interval [—0.01,0.01]. In the 
early stages of the fracture process, the stress field is cal- 
culated in the presence of these variations until buckling 
commences. Before sizable cracks appear, forces combine 
to flatten out the vertical displacements. That is, a flat 
configuration is energetically preferred to begin with, and 
fracture propagates according to the non-buckling simu- 
lation. As significant cracks begin to appear, however, 



the conditions at some point become favourable for the 
out-of-plane components of the stress field to be real- 
ized, and buckling sets in. From here on the random 
noise is discarded, and the next displacement configura- 
tion is simply calculated from the previous coordinates. 
When the lattice has been broken, a new set of vertical 
displacements is generated for the next sample, i.e., a 
sample-specific random noise is used. 

Lattice buckling modes in the presence of a center- 
crack of size ~ L/3 are shown in Fig. ^0 The only dis- 
order present here is that due to the out-of-plane initial- 
ization, but evidently a number of buckling modes may 
appear. In the following, cases where the deflection of the 
edges of the crack is to the same side, i.e., up-up or down- 
down, are referred to as symmetric buckling, and cases 
where the deflection is to opposite sides is referred to as 
anti-symmetric buckling. In Fig. 1131 (D) and (E) are ex- 
amples of the former, and (F) and (H) are examples of 
the latter. Another buckling mode which frequently ap- 
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FIG. 14: A beam lattice of size L — 50, showing (a) the 
symmetric and (b) the anti-symmetric buckling modes for a 
crack between I = 12 and / = 40. 



pears is that shown in (A), (G) or (J), where the main 
bulge at one of the crack edges is made up of four, rather 
than three, half-waves. Based on the few dozen samples 
observed, (A) and (J) evolve into type (D) after a few 
more beams have been broken, and (G) into type (E), 
i.e., the symmetric buckling mode prevails in each case. 
However, due to the randomness introduced, examples 
such as (H) are not completely anti-symmetric about the 
neutral plane. Hence, even for this simple crack config- 
uration, the exact shape of the out-of-plane deflection 
can vary considerably. The overall shape, however, tends 
to fall within the main categories, i.e., one of the two 
symmetric or anti-symmetric buckling modes. This way 
of initializing the out-of-plane deflection is suitable for 
studying disordered systems, where a lattice without any 
initial geometrical discontinuity is strained until random 
cracks begin to appear. As the cracks grow buckling sets 
in at some point, depending on the configuration, posi- 
tion and size of the initial cracks. 

Another way of initializing the out-of-plane deflection 
is to impose a small vertical deflection on a very few nodes 
in strategic positions. This is most practical when study- 
ing an "ideal" buckling scenario, such as the fracture of 
a non-disordered plate with a perfect center-crack. 



VII. FRACTURE CRITERION 

In order to study how buckling affects the fracture 
properties of a two-dimensional structure an appropri- 
ate breaking criterion should be chosen. As previously 
mentioned, this can be done to suit a range of engineer- 
ing requirements. Often the mode of rupture in the out- 
of-plane direction is radically different from that which 
takes place within the plane. In paper or cloth, for in- 
stance, the phenomenon which first springs to mind is 
tearing. The energy required to propagate a crack across 
a given area in tear mode is much less than that which 
causes the same area to fracture in pure tensile loading. 
This is especially the case with paper. 



Out-of-plane contributions to the breaking criterion 
must be included by some other mechanism than that 
provided by Eq. ©, since the latter is relevant to re- 
gions which are comparable in size to a beam. The stress 
intensification due to buckling, on the other hand, is due 
to much smaller regions, i.e., comparable in extent to the 
sharp crack tip. One way of enhancing the stress due 
to buckling is to combine torsion with axial stress. The 
larger the load, the more sensitive the beam will be to 
the presence of a given amount of torque. Compressive 
loads are assumed to alleviate the torsional moment, but 
only to a very small degree. 

Hence, the breaking criterion can be stated as 
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is the enhancement factor in Eq. I)42[l. 

Considering Fig. 1141 the breaking stress is increased in 
case (b) and also in case (a) provided the deflections in 
front of and behind the crack are not congruent. When 
the bulges are completely symmetric, however, (a) does 
not intensify the breaking stress and thus contradicts ex- 
perimental findings. 

Another possibility is to assume a crack-tip stress en- 
hancement which depends on the out-of-plane bending 
moment. Here the in-plane displacement component yi 
observed in Fig. [5J i.e., the backward and forward move- 
ment of the crack edge, creates an angular displacement 
about the X-axis. For a sufficiently thin plate the resis- 
tance towards bending will not be sufficient to halt the 
out-of-plane deflection once it has commenced, since the 
forces involved act over a region much larger than the 
immediate neighbourhood of the crack. Due to the short 
distance which separates the top and bottom surfaces, 
the resulting "lever-arm" effect creates an asymmetric 
stress-gradient across the crack front in the direction of 
the thickness. Whereas tensile force on the concave side 
is then reduced, it increases on the convex side. This 
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increase comes in addition to the stress already concen- 
trated along the crack front, i.e., the very presence of a 
crack creates a screening effect which re-distributes the 
in-plane stresses so as to cause a build-up in the load 
at the crack tips. For a crack that has grown to an ex- 
tent which allows buckling to occur, this in-plane stress is 
significant. The crack-tip opening angle also plays an im- 
portant role. Buckling in brittle materials, for instance, 
is known to have a profound effect on the maximum load 
the system can tolerate before breaking. In a FEM study 
by Seshadri and Newman [T^ a hypothetical very large 
critical crack-tip opening angle was used to model buck- 
ling in a ductile material. Strength reduction in this case 
was found to be significantly smaller than for brittle ma- 
terials. 

In the beam model, the crack tip is never sharper than 
exactly one beam length. To emulate the above stress en- 
hancement due to out-of-plane bending we instead imag- 
ine a sharp crack to be embedded within that beam which 
on the lattice defines the tip of the crack, and consider a 
combination of axial stress and moment. Out-of-plane 
bending modes are shown in Figs. El an d El where 
the displacements of the schematic lattice at the top 
have been exaggerated somewhat to illustrate the point. 
Specifically, the Z-displacements of contour B have been 
scaled up 100 times with respect to those of contour A, 
which itself is scaled up with respect to the horizontal 
extent of the lattice. The in-plane ^-displacements have 
also been adjusted accordingly. 

Experimental evidence indicates that the stress en- 
hancement at the crack tips is more or less similar in 
the symmetric and anti-symmetric buckling modes. To 
incorporate this we distinguish between the two cases. 
Hence, retaining Eq. (|41|1 . we introduce 
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FIG. 15: Symmetric buckling. Shown at top, for a lattice 
of size L = 100, is (A) the / = 51 contour, passing through 
the middle of the lattice, and (B) the I — 34 contour, pass- 
ing through the left-hand side crack tip. Also shown is the 
bending mode of the beam which defines the crack tip at the 
junction between I — 34 and J = 51. At the bottom are 
shown the out-of-plane angular displacements Vi about the 
X-axis, in the case of (B) above. Also shown are the neigh- 
bouring contours, (C) I = 33 and (D) I = 35, where (D) is 
discontinuous due to the intersecting crack. 



and 



F c = - max( 
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(51) 
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replaces in Eq. (|42|) . In this prescription, 



M. 



U) 



±M. 



(48) 



denotes symmetric (— ) or anti-symmetric (+) buckling, 
respectively, with the signs referring to the direction of 
the moment at the two beam ends. 

For the effective stress in the beam, we now use 



= - X 



MP - Mf 
t j 



in the symmetric case, i.e., when 



M 



W 



-M, 



(i) 



(49) 



(50) 



in the anti-symmetric case, i.e., when 

1 J 

The enhancement factor in Eqs. I|49fl and (|51fl is 



A; 



X 

where 



' 2 L 
0, 



U) 



L 



< 0, 

fP > o, 



A l = piA i:X + qiA. 



(52) 



(53) 



(54) 



is a discontinuity operator. The choice made above 
causes the breaking stress of the beams at the crack tips 
to increase by a comparable amount in symmetric and 
anti-symmetric buckling. 

The expressions for x an d X m Eqs. I)45|) and 
respectively, have been chosen, very generally, to incor- 
porate some overall effects related to size, material and 
relative dimensions. Hence, it is reasonable to assume 
that, for a given size, "tearability" , or the "lever-arm" 
effect, increases with decreasing sheet, or plate, thick- 
ness. In conjunction with this, the crack-tip opening an- 
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FIG. 16: Anti-symmetric buckling. Shown at top, for a lattice 
of size L = 100, is (A) the / = 51 contour, passing through 
the middle of the lattice, and (B) the / = 34 contour, passing 
through the left-hand side crack tip. Also shown schemati- 
cally is the the beam which defines the crack tip at the junc- 
tion between / = 34 and J = 51. Although the bending 
mode is correct, the actual angles at the ends of the beam are 
negative and not positive as shown. Out-of-plane angular dis- 
placements Vi about the X-axis are included at the bottom, 
with the notation being the same as in Fig. 1151 



gle, which decreases with increasing resistance towards 
in-plane bending, also enters the picture. Both effects are 
presently included via the ratio of the in-plane to the out- 
of-plane inertial moment for bending, i.e., a 2 — Iz/Ix- 
The length of the arm with which the out-of-plane forces 
act is assumed to be proportional to the vertical extent 
of the buckling zone. Since this, in turn, is proportional 
to system size, a factor L\f£ | is also included. As noted 
previously, fracture is initiated by displacing the top row 
of the lattice a fixed distance, usually corresponding to 
one beam length. To avoid scale effects associated with 
this, a further factor Lq is included, where Lq is the size 
of the reference system for which the top row displace- 
ment is exactly one beam length. The introduction of 
a reference system allows for the possibility of compar- 
ing systems of varying size where the physical behaviour 
involved requires the same relative external boundary 
conditions. For instance, referring to the intact lattice, 
mode-I loading then imposes the same initial strain of 
(Lq + 1) _1 on each beam. Although computational time 
increases when L > Lq, features such as how various 
buckling modes appear with respect to system size will 
depend on the external loading. Otherwise, Lq = L 
might probably be used in cases where we are interested 
in features which depend on the internal processes of the 
fracture mechanism, such as the roughness exponent of 
crack interfaces. 



To guard against unphysical breaks, we introduce 

PyA — n y : i-l + (55) 

which contributes when one or both nearest lateral neigh- 
bours are intact, and 



c L , m -i 

3=1 
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)+ II ^""iM+j) ( 56 ) 

3=1 



which contributes when a certain number of neighbours 
have been broken. For any node i, the array 
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(57) 



now keeps track of the status of the beam which extends 
away from i in the direction of the y-axis, i.e., it remem- 
bers whether this is broken or intact, respectively. The 
combined expression, 



has the property 



— ^~*yA ' ^3y,*> 



^■y,i - -I i 



(58) 



(59) 



as has K x _i governing cracks in the normal direction. In 
other words, Eq. (|54|1 ensures that the stress enhance- 
ment mechanism is activated only in cases where the lat- 
eral neighbour on one side is intact while simultaneously 
a certain number of beams, defining a minimum crack 
length CL, m , are broken on the other side. 

In most cases the operator Aj is not necessary. It 
has been included to avoid cracking being induced near 
the top and bottom rows of the lattice. For very large 
systems, and especially in cases where L is significantly 
larger than Lq, breaks sometimes occur due to the large 
angular gradients in beams extending up from J = 1 or 
down from J = L+2, see, for instance, Figs.ll5landllb1 In 
the present formalism the properties of beams with inflec- 
tion points are not considered, the smallest crack that can 
cause buckling, i.e., a bulge consisting of at least three 
half-waves, is therefore approximately C^.m = 4. This 
is confirmed in numerical runs for systems with small L, 
but where L >> Lq. A problem with using such a large 
value of Ch,m is that it excludes cracks inclined at an 
angle with respect to the horizontal. Over a wide range 
of system parameters and external boundary conditions, 
however, Cl,™, = 2 was found to be adequate. 

In the limit of no buckling, i.e., Sv — > or Su — > 0, 
Eq. (|41|l reduces to Eq. ©. In other words, if buckling 
is not activated, then neither is the stress enhancement 
mechanism in the fracture criterion. 

Finally, in order to illustrate how Eq. 141fl . with 
Eqs. (|46|l to l|59|) defining the stress enhancement mech- 
anism, works within our model of buckling, we consider 
the buckling response ratio of the residual strength of the 
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FIG. 17: The buckling response ratio, Ft, shown as a func- 
tion of the crack-length-to-thickness ratio, Ct/t, for various 
systems with L — 3Cl- Open squares denote the results of 
varying L while keeping t constant, filled circles denote the 
results of varying t with L = 54 fixed, and crosses denote 
similar results with L — 24 fixed. Kuhn and Figge's linear 
expression [l^| is also included for comparison. 

system. That is, F\, = A /A^, where Ao and A z repre- 
sent the maximum applied external force a restrained or 
buckled plate, respectively, can tolerate before breaking 
apart. Early experimental results show that the decrease 
in strength due to buckling increases as the ratio of the 



crack-length Cr, to the thickness t is increased. A linear 
relationship was proposed by Kuhn and Figge which, 
in the case of brittle materials, has been shown to agree 
well with more recent FEM calculations [13 • In Fig.EI 
results obtained with the beam model are compared with 
the Kuhn-Figge relationship. A small correction to the 
size of the central crack has been made to account for 
the finite size of the beams. The effect is very small, 
shifting the values of the smallest systems slightly to the 
left, thus improving the agreement with the Kuhn-Figge 
relationship from very good to excellent. 



VIII. SUMMARY 

To summarize, we have included the additional degrees 
of freedom necessary to describe the interaction of cracks 
with buckling in the elastic beam model. This model is 
stochastic in nature, so that sheets with random cracking 
at any level of meso-structural disorder can be studied, 
including systems with no disorder. In addition to impor- 
tant issues of practical relevance in traditional fracture 
mechanics, such as strength properties and stability, the 
present model also enables fundamental aspects of frac- 
ture in random media to be explored. 
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